The energy-scale-dependent Composite Operator Method 
for the single-impurity Anderson model 



o 
o 



oo 



o 
o 



> 

(N 
00 

in 
o 

O 

-i— > 



c 

o 
o 



X 



Adolfo AvelljQ and Ferdinando Mancini 
Dipartimento di Fisica "E.R. Caianiello" - Unita di Ricerca INFM di Salerno 
Universitd degli Studi di Salerno, 1-84081 Baronissi (SA), Italy 

Roland Hayn 

Institut fur Festkorper- und Werkstofforschung (IFS) Dresden, D-01171 Dresden, Germany 

(Dated: February 6, 2008) 

The recently developed energy-scale-dependent Composite Operator Method is applied to the 
single-impurity Anderson model. A fully self-consistent solution is given and analyzed. At very low 
temperatures, the density of states presents, on the top of the high-energy background, a Kondo- 
like peak whose parameter dependence is discussed in detail. The proposed method reproduces the 
exact results known in the literature with very low numerical effort and it is applicable for arbitrary 
values of the external parameters. 
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I. INTRODUCTION 

The theoretical description of strongly correlated elec- 
tron systems, like transition-metal oxides 1 - and heavy 
fermion compounds 2 , is of high actual interest and ef- 
fective analytical methods to study them are looked for. 
A general problem consists in connecting the high- and 
low- energy scales in a proper way. Methods based on the 
use of the equations of motion for the Green's functions 
(e.g., the projection methodo 3 : 4 : 5 : 6 : 7 : 8 : 9 : 10 : 11 : 12 : 13 ), usu- 
ally give a rather reliable description of the high-energy 
features of strongly correlated systems, but, quite often, 
the overall solution obtained by their application does 
not reproduce the low-energy physics accurately enough. 
Other techniques, like the slave-boson approximation 14 , 
provide a correct picture at low energies but fail at 
higher ones. However, there is an emerging consensus^ 
that in strongly correlated electronic systems we should 
search for a correct description of both energy scales at 
once. The anomalous behaviors shown by these systems 
are caused and/or influenced by both the broad inco- 
herent spectral features far away from the Fermi level 
and the more dispersive quasi-particle bands close to it. 
Recently^, it was shown how the Composite Operator 
Metho d 12 ' 13 (COM) could be used to resolve coherent 
low-energy features in a proper way by solving the SU(N) 
Kondo model. 

The Kondo model should be considered just as a first 
step. It contains no charge degrees of freedom of the 
strongly correlated impurity level, and it is usually not 
sufficient for a realistic description of 4/-electron spec- 
tra. For this latter purpose the single-impurity An- 
derson model is much better suited^. Furthermore, 
the solution of the Anderson model is a building block 
within the Dynamical Mean Field Theory^ (DMFT) al- 
gorithm, which has favored large progresses in the com- 
prehension of the Mott metal-insulator transition phe- 
nomenon. The single-impurity Anderson model is well 



and widely studied^ 7 -, but the reliable methods are rather 
involved (e.g., quantum Monte Carlo method™, Bethe 
ansatz 19 , numerical renormalization groupSi and non- 
crossing approximation^ 7 ") and require a quite huge com- 
putational effort. This is certainly an obstacle to use 
them within the DMFT algorithm or to interpret the re- 
alistic spectra of rare-earth compounds. According to 
this, we present here a simple analytical method to solve 
the single-impurity Anderson model, which is capable to 
reproduce both the high- and the low- energy features, 
as known after the exact solutionis, in a reasonable way 
and, practically, without requiring any computational ef- 
fort at all. 



II. HAMILTONIAN 
AND EQUATIONS OF MOTION 

The single-impurity Anderson model is defined by the 
Hamiltonian 



k.tr 



V 



k.cr 



where c^ a represents the electrons in the valence band 
(ek) and f a those at the impurity level (e/); rif a is the 
density charge operator for the /-electrons of spin a, N 
is the number of sites, U is the Coulomb repulsion at the 
impurity level and V is the strength of the hybridization 
between the valence band and the impurity level. For the 
sake of simplicity, we restrict ourselves to a simple two 
(spin) degenerate impurity level (er = ±1, a = — a). 

The first step within any method based on the projec- 
tion technique consists in individuating an appropriate 
set of basis operators. In fact, for strongly correlated sys- 
tems the original electronic operators are just the wrong 
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place where to start any approximate treatment 1 . The 
electrons completely lose their identities owing to the 
strong interactions and other complex excitations appear. 
These latter are the only effective quasi-particles present 
in the systems and only in terms of them any description 
of the dynamics should be attempted. The basis (com- 
posite) operators should be then chosen in order to de- 
scribe such excitations. The choice is not easy, but some 
recipes could be givenS^. In the present case one should 
have as minimal requirement an appropriate description 
of the free propagation of c-electrons and the atomic dy- 
namics of /-electrons. The latter have two possible exci- 
tations corresponding to the transitions from |0)/ to |cr)/ 
and from to | TJ.)/ ■ That leads to the first three basis 
operators 



01<T = C 0a 



-y 



02cr = £<r = (1 - nfff)f a 
03<r =r\a= nfsja . 



(2.2) 



It is worth noting that 02<r — £,a = X Ua and 
aX a2 are the Hubbard operators describing the lower- 
and the upper- Hubbard subbands. 

These three operators give a very good description of 
the high-energy features, but are certainly insufficient to 
describe the low-energy ones and the mixing between the 
two regimes. Then, in order to overcome such limitations, 
we write down the equations of motion for the operators 
(l2~2l 



basic excitations. Such basis cannot be considered com- 
plete according to the fact that an infinite degree of free- 
dom system possesses an infinite number of operator ba- 
sis. Then, it is obvious that the basis is neither unique. 
In particular, it has been already mentioned above that 
many recipes can be used in order to construct an oper- 
atorial basis, each according to the features of interest. 
However, as we will see from the obtained results, this 
basis seems sufficient to catch the main features of the 
dynamics of this system. 

In order to be more compact, we rewrite the six basis 
operators in spinorial notation 



02 = £ = (1 - n f )f 

03 =T) = Tiff 



04 = |(1 - n f )c 

05 = \a on s ■ c 

06 = 4 • £ ® v ■ 



(2.6) 



Here ftt = /t • a • f, a are the Pauli matrices, • denotes 
the scalar product in spin space, o the scalar product in 
direct space and (g> the tensor product. This notation 
will be used hereafter. It should be noted that the ex- 
act expressions of the basis operators have been chosen 
such that they transform under the particle-hole trans- 
formation (co — > cj and / — » /t) into themselves or into 
another basis operator: 0i — > ip\, 02 — ► 0|, 03 — * tpl 
and 0„ — ► — 0J, for n = 4, 5, 6. 

As next step, we derive the equations of motion for the 
retarded Green's functions (GF) 



V 

idtCka = [Oka, H] = e k Ck CT + —7=(€o- + V<t) 

v N 

V 

v 

id t r] a = [r)a,H\ = [e f + U)rj a + —c 0a - Vn a 



(2.3) 



and account for the appearance of a new operator, 
namely the fluctuation field 



04o- + 05cr + 0. 



6<T 



(2.4) 



with 



04<T = ^(1 - rif )c cr 

05O- = crc QcT S z + c 09 S° 

06ct = c ogfsf<T t 



(2.5) 



where n/ = n/ T +n/;, S z = (n/ T -n/;)/2 and^ = /t/ ff . 
The fluctuation field tt describes the coupling of the va- 
lence band to density-, spin- and pair- impurity fluctua- 
tions and opens the possibility to study the low-energy 
dynamics connected with them. According to this, we 
have decided to include also the operators l|2.5|) into the 
basis and study the system in terms of this set of six 



= T[6(t-t')({Mt),^ m (t')})] 



(2.7) 



where T stays for the Fourier transform and TZ for the re- 
tarded time ordering operator. The inhomogeneous term 
in the equations of motion for the Green's functions G nm 
is constituted by the normalization matrix (its expression 
in terms of basic field correlators is given in Eq. (|A1|) of 
App. A) 



({0„(i),0t l (t)}) , 



(2.8) 



where { — ,— } stays for the anticommutator and (— ) for 
the thermal average. The normalization matrix is not 
only fundamental to derive the Green's functions but pro- 
vides also important information about the total spectral 
weights of the fields. Its matrix elements depend on the 
expectation values Cy = (0i0>J) (correlation matrix). It 
contains for instance the average charge density at the 
impurity level (n/) = 2(1 — C22 — C33), and the dou- 
ble occupancy Df = (n/fn./j) = I — C22 — 2C33. For 
further use we define also the matrix elements with the 
fluctuation field, i.e., I nw = ({ip n (t) , (t)}) . Then, af- 
ter the equations of motion l|2.3|l . we have the following 
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expressions for the Green's functions 

Gn - r + V 2 T 2 (G 22 + 2G 23 + G 33 ) 

B+(C+ + x 



G 



22 



G 



£4 

~F 



I33- 



G23 — X-G22 — 



F 2 
B 2 _ 

VT 



Gjj 



-G_ 



Gtj 



F 2 

B_(G_ +x_G 4 



F 2 

( J B+G 7r7r + G++a: + G_) 



(2.9) 



where we introduced the following abbreviations 



1 



Xs = ±v 2 r r s 

c+ = (J 2w + ±vt /i.)i% 

F = 1 — x + a;_ 



r — 1 

- 1 - s LO-e 3 -V*T /2 

B s = VT S - Vx s Y s 

C- = (I 3w + ^VT hn)T- 

(2.10) 



with s = ±, s = — s, e+ = Ef, and e_ = Sf + U, to 
shorten the notation. It is worth noticing that Tq(uj) is 
just the free (V = 0) propagator for Co operator. Wc 
see from (12. 911 that we can calculate the Green's function 
for the valence electrons Gn and the /-impurity G// = 
G22 + 2G23 + G33 once we know the Green's function of 
the fluctuation field G n7r . 



III. HIGH- AND LOW- ENERGY SCALES 



with 



(3.4) 



After Eq. 13.11 and the above reported consequent 
reasoning, we split the still unknown Green's function 
Gtttt = T(K [7r(£)7r(i')]) into a low-energy component 
G^ = T{U [n L (t)ir L (t')]) = F{U [^{t)^{t')}) = Gg 
and in a high-energy one G^ = T(K [tt h (t)?: 11 (t')] ) = 
EJU=4^(ft [1>S (*)C(0]>- Wc neglect the cross-term 
T(TZ [tt l (t)ir H (t')] ) according to our previous assump- 
tion of no-overlap of the energy sectors spanned by if>^ 
and all other fields involved. According to Eq. (|3.2|) , Gg 5 
obeys the following equation of motion 



LOG 



55 



J 55 



KlGi 5 



K 3 G 35 



(3.5) 



where Gf 5 = J 7 {TL [^(t) ^5 (*')])• Then > for the sake 01 
consistency, we approximate Gf 5 (see appendix for com- 
plete expressions) only by those components that are ex- 
plicitly proportional to G§ 5 , i.e., 



@ib - ^ r o(G2 5 + G% 5 ) 

^25 — "^^55 
^35 — —pT 1 * 55 



(3.6) 



In order to resolve the low-energy features embedded 
in the high-energy background, following the idea given 
in He:'. Hi, which is based on the well-established physical 
assumption that at low energies we have a quasi-particle 
theoryii as also derived by the slave-boson theory^, we 
split the dynamics of the fluctuation field into an high- 
and a low- energy part. As the essence of the Kondo 
effect consists in the coupling of the valence band to the 
spin fluctuations at the impurity level, we split only ip5 



05 = *l>6 + 1>> 



(3.1) 



and assume instead that the charge (^4) and the pair (i(jq) 
terms in the fluctuation field are sufficiently well repre- 
sented by their high-energy parts only (04 = ip^ and 
06 = 06 )• I n practice, we assume that the low-energy 
field i/jL spans a different energy sector of the Hilbcrt 
space with respect to 0f0 04, 06, an d describes a coher- 
ent quasi-particle at very low energies; energies that are 
much smaller than any other defined in the Hamiltonian. 
According to this, we make the following ansatz 



idtipi = [05, #] = Kicd- 



+ K3V 



(3.2) 



The coefficients (i — 1,2,3) are determined by pro- 
jecting onto the basis 12. 6|) 



«i 



V(L 



I35) 



«a = efl&fln + VI^/{2I 22 ) 

k 3 = (e f + U)l£ 5 /I 33 -Vl£ 5 /(2I 33 ) 



(3.3) 



and we finally obtain 



^TTTT _ ^55 ~~ 



T L 



On 



(«iVT 



■-n 

K 2 )S 4 



(kiVT + k 3 )b. 



(3.7) 



2F 



It is worth noticing that, in determining the actual ex- 
pressions for Gf 5 , we have neglected both the cross-terms 
F{K [ip?(t)ip£ (*')]) with i = 4,5,6, according to our 
previous assumption of no-overlap of the energy sectors, 
and those terms proportional to T and r s , according to 
the well-defined high-energy character of these latter. 

We are now left with the task of computing the 
high-energy contribution G%„ = El!, m =4 G L where 

Gnrn = F(R> [V'n (*)^m (*')]) • In order to accomplish 

this task, we use the mode-coupling approximation 23 , 
also known as self-consistent Born approximation. In 
practice, we neglect the mixing terms among different 
bosonic modes (i.e., we take Gf 5 — G|g — Gf| = 
0) and decouple the remaining corresponding time- 
ordered propagators in terms of the charge (So = 
_)), spin (S z = F{T[S z {t)S z {t')])) and 

fl{t)fr{t)f\{t')f\ (*')])) impurity time- 
ordered propagators, and of the valence electron time- 
ordered propagator (S11). For instance, we have 



T{T[n f (t)n f {t r 
pair [S p = J-(T 



SIM = 3^ 



dams»(w-n)5ii(n) 



(3.8) 
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Actually, for energies as high as those we intend to de- 
scribe in the high-energy sector and as far as the impu- 
rity level is not too deep inside the valence band (i.e., 
for \ej\ not too small with respect to the bandwidth) we 
can safely take the atomic limit for the spin, charge and 
pair impurity propagators. In the high-energy regime, 
any other treatment would not substantially affect our 
results (i.e., by taking the atomic limit we already get 
an accuracy of the same order of the energy scale we are 
computing). For instance, we have 



S s (u) = --2m[(n f )-2D f ]6(w) 



(3.9) 



Then, in the symmetric case (— 2er — U, leading to 
(nf) = 1), all the non-zero contributions are proportional 
to Gn, namely 



G»=±(l-(n f )+2D f )G n 
Gg = 3 -((n f )/2-D f ))Gn 

G 66 = -DfGn 



(3.10) 



Finally, we can write the following expression for the 
fluctuation field Green's function G TT , which take into 
account both high- and low- energy contributions 



G>xiX 



2^)611 + . 



(3.11) 



For Ef outside the valence band we have numerically 
checked that the influence of the double occupancy D f in 
Eq. i|3.11[) is very small and could be neglected. This lat- 
ter procedure has been adopted in the asymmetric case 
(as U has always been taken much larger than the band- 
width) with the result that Eq. (|3.11ll . with Df = 0, also 
holds. 

It is worth noting that the high-energy Green's func- 
tion G^ 5 exhausts only part of the spectral weight J 55 , as 
expected from an independent evaluation of normaliza- 
tion matrix, 



l imG * 3(^£> 



lim G 55 



(3.12) 



which can be derived from l|3.10|l using lim w> i G\\ = 
\/u). On the contrary, we have: lim w> i wG^ = la — 
lim w> i G 44 and lim w> i loGqq = I 66 = lim a)> i G 66 . This 
consideration is fundamental as explains why we are al- 
lowed/forced to consider the possibility of a low-energy 
contribution to the dynamics and why we can/have to 
concentrate on the spin fluctuations. 



IV. SELF-CONSISTENCY 

In order to finally compute the GF in a self-consistent 
way, one has to take into account two aspects. First of 



all, we have to determine the matrix elements of the nor- 
malization matrix by calculating the expectation values 



duj[l - / F (w)]Ay(cj) 



(4.1) 



where Aij(uj) = — i3[G^(w)] and f-p is the Fermi func- 
tion. Second, a self-consistent calculation of the low- 
energy parts of the spectral weights , I35 , and I$ 5 is 
needed in order to determine Gg 5 (see Eq. (13.611 . For that 
purpose we start from the total spectral weights I25 — 
-3Gi2/2, I35 = -3Gi 3 /2, and I 55 = \ {{n s )-2D f ) + 
G15. Then, for a given temperature T and for given 
low-energy contributions I L , one can determine the to- 
tal weights /, written schematically as I[T,I L ), by ex- 
ploiting the well-known connection between the correla- 
tion and the Green's function matrices l|4.1ll . To extract 
the weight connected with the low-energy part, we com- 
pare the spectral weights calculated for I L = (defin- 
ing I H [T] = I[T,I L = 0]) with the complete expression. 
Then, the required self-consistency condition can be writ- 
ten as 



I L [T] = I[T,I L }~ I H [T] 



(4.2) 



This procedure has demonstrated to be capable to grasp 
part of the effects connected with the mixing between 
the high- and low- energy sectors. To calculate G12, C13, 
and G15 we need also the corresponding Green's func- 
tions whose expressions are given in the appendix for 
completeness. 



V. NUMERICAL RESULTS 

The numerical calculations have been performed with 
a density of states for the valence band electrons c(k) 
which is constant and finite only between energies —D 
and D. That is, we have used a free (V = 0) propagator 
To (a;) for cq with the following well-known expression^ 



r (w) 



1 

2D 



In 



D - fj, 



2D 



0(D-\uj + n\) (5.1) 



We have given the formula for an arbitrary value of the 
chemical potential fj,. The expectation values at finite 
temperature T have been calculated by means of the Mat- 
subara formalism 



C lm = -2TM[J2(G 



hn {^n ) 



n=0 



:)] + [!- Me)]I tm 
(5.2) 



where uj n = (2n + \)irT . The arbitrary small number < 
e«1 and the Fermi function /p have been introduced in 
order to ensure the correct behavior at high frequencies. 
The zero-order Green's function along the imaginary axis 



5 



is given by 



1 (D + rf+u* 



N,(co) 



Z / U) 

+ ~ arctan 



2L> 



D + n 



arctan ■ 



D - fi 



(5.3) 



Let us start the discussion regarding our results with a 
short briefing about the high-energy features. It is worth 
noticing that our solution manage to describe the physics 
at the energy scale of U in a very effective way by explic- 
itly taking into account the excitations related to the 
Hubbard subbands. In Fig. ^ m fact, we can clearly see 
that the relevant features are correctly reproduced: the 
splitting of the non-interacting band into two subbands 
for finite values of U, a distance in energy between the 
centers of mass of the two subbands almost equal to U, 
spectral weights independent of bare U, but strongly de- 
pendent on the ratio V/U. As regards the scale of energy 
V, we use a basis rich enough to give the exact solution 
in absence of U : second order (V 2 ) resonant behavior, 
spectral weights satisfying the ordinary sum rules. Then, 
higher order processes ruled by the competition of these 
two scales of energy (the localizing U and the dispersing 
V), and strictly dependent only on the ratio V/U, are 
just those responsible for the low-energy features we will 
discuss in the following paragraphs. 

We can now move to the study of the low-energy 
features (whose description is the actual goal of this 
manuscript) and, in particular, to the analysis of the evo- 
lution of the Kondo peak by moving the impurity level to- 
wards the valence band in the symmetric case (see Fig.[TJ. 
The peak is widened for decreasing values of |e/| outside 
of the band. In that case we have two bound-states at 
roughly ±Ef which are not shown in Fig. ^ For \sf\ in- 
side of the band we can observe a characteristic three 
peak structure. Remarkably, the central peak height is 
not changed and identical to the height for U = 0, in 
agreement with the Fridel sum rule. All these features 
are in agreement with the exact behavior known from the 
numerical renormalization groupit 

Fig. 2 shows results in the asymmetric case. The 
Kondo peak is shifted slightly above the Fermi level and 
its shape becomes asymmetric. Furthermore, the width 
of the Kondo peak is decreased in agreement with the 
reduction of the Kondo temperature. It is worth not- 
ing that the Kondo peak also remains for U — > oo. If 
we increase the temperature, the Kondo peak vanishes 
by diminishing its width (not shown). We have de- 
fined the Kondo temperature as the one at which the 
parameter I§ 5 has a change in the concavity when plot- 
ted as a function of the temperature. The Kondo tem- 
perature becomes exponentially small for small values 
of J = V^ 2 (l/|e/| + l/\ef + U\) in agreement (see Fig. 
3) with the formulail T K = D^/2jp^exp[-l/(2Jp )], 
where po — l/(2D) is the density of states of the un- 
perturbed valence band at the Fermi level. It is worth 
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FIG. 2: Kondo peak in the asymmetric case with the param- 
eters e s = -2, V = 0.7, D = 1, fi = 0, and T = 0.001. The 
width of the Kondo peak diminishes with increasing Hubbard 
correlation U [lt 5 = 0.374 for U = 4 (symmetric case), 0.287 
for (7 = 6, 0.215 for U = 8 and 0.068 for U = col. 
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FIG. 3: Kondo-temperature in the symmetric case in depen- 
dence on the exchange coupling J compared with the formula 
given in the text. 
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noting that this formula is correct only in the symmet- 
ric case and for J -C D. We have a strong dependence 
on the value of the hybridization V that shows a lack 
of universality. Anyway, the linear behavior at small J 
is preserved for any value of V showing that the expo- 
nential dependence, which is not possible to obtain per- 
turbatively, is correctly described. We also studied the 
infinite bandwidth case, but it does not give any quali- 
tative change to the physical picture described above. 

In conclusion, a recently developed energy-scale- 
dependent approach**, which is capable to reproduce in 
a reasonable way both high- and low- energy features 
of known exact solutions of impurity models, has been 
extended to the case where relevant charge fluctuations 
are present. The originating procedure, the Composite 
Operator Method****, provides a fully self-consistent so- 
lution where, on the top of a broad high-energy back- 
ground, a Kondo peak is present at low temperatures. 



The parameter dependencies of the peak features and of 
the Kondo temperature have been correctly reproduced, 
with respect to the exact results known in the literature, 
with very low numerical effort. 
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APPENDIX A: THE NORMALIZATION MATRIX 

The normalization matrix is found to be 
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APPENDIX B: THE SELF-CONSISTENCY where we have in close analogy to (jo~TT)> 

CYCLE 

To close the self-consistency cycle, the following GF 
are needed 



G12 — VTq(G22 + G23) G13 — VTq(G23 + G33) , 



3 f(n f ) 



(Bl) = - I i-LL - D f ) G u + G^ 5 (B3) 



and 



Gib - VT (G 25 + G 35 ) 

G 25 = ^G. 5 + (C5++ ; +C5 - } (B2) 

G35 = -—^rG^ + + ^~ C5+ ^ and G 5+ = / 2 5r+, and G 5 _ = / 35 r_. 
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